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Abstract. We present a new multi-fluid, grid MHD code PIERNIK, 
which is based on the Relaxing TVD scheme. The original scheme has 
been extended by an addition of dynamically independent, but inter- 
acting fluids: dust and a diffusive cosmic ray gas, described within the 
fluid approximation, with an option to add other fluids in an easy way. 
The code has been equipped with shearing-box boundary conditions, 
and a selfgravity module, Ohmic resistivity module, as well as other 
facilities which are useful in astrophysical fluid-dynamical simulations. 
The code is parallelized by means of the MPI library. In this paper 
we shortly introduce basic elements of the Relaxing TVD MHD algo- 
rithm, following Trac & Pen pUDg)) and Pen et al. J20S31), and then 
focus on the conservative implementation of the shearing box model, 
constructed with the aid of the Masset's (|2000j) method. We present 
results of a test example of a formation of a gravitationally bounded 
object (planet) in a self-gravitating and differentially rotating fluid. 

1 Introduction — the basic Relaxing TVD scheme 

The Relaxing— TVD conservative scheme (Jin & Xin ,1995;) presented by 
Trac & Pen (|2003p and Pen et al. (|2003p . who provided short codes with a basic 
implementation of the method, is a second order algorithm in space and time. 
The scheme efficiently deals with shocks without artificial viscosity. The code is 
very flexible and can be extended with new modules representing additional phys- 
ical processes. The simplicity and robustness of the code is reflected in general 
performance of 10 5 zone-cycles/s (on 2 GHz AMD Opteron processors). 
The conservative form of MHD equations serves as a starting point 

d t u + d x F(u, B) + d y G(u, B) + 9 z H(u, B) = 0. (1.1) 
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where u = (p 7 m x ,m y ,m z , e) and F(u, B), G(u, B), H(u, B) are fluxes of the 
conserved fluid variables in x, y and z directions, respectively. Other symbols 
have their common meaning i.e. p denotes density of the fluid, whereas m x , m y , 
m z are components of momentum density and e stands for total energy density. 

Firstly, a dimensional splitting of equation (jl.ip is made by constructing nu- 
merical solution with timestep At to the equation 

d t u + d x F(u, B) = 0, (1.2) 

separately for each dimension. Then u and F are split into waves that propagate 
leftwards and rightwards 

= \ (u - |Y u* = \ (u + \ \ , u = + u* (1.3) 

F L = cu L , F R = cu R , F = F R + F L (1.4) 

where c, called the freezing speed, is a function that satisfies c > max(|i> ± Cf\), 
v is the fluid speed and c/ is the fast magnetosonic speed. Now, equation (|1.2|) is 
equivalent to two independent equations: 

d t u L - d x F L = 0, d t u R + d x F R = 0. (1.5) 



The above pair of equations ()1.5j) is solved by means of an upwind conserva- 
tive scheme, separately for right- and left-going waves, using cell-centered fluxes. 
To achieve 2nd order spatial accuracy, a monotone upwind interpolation of fluxes 
onto cell boundaries is made, with the aid of a flux limiter. Second order accu- 
racy of time integration is achieved through the Runge-Kutta scheme (details see 
Trac & Pen [20031 Pen et al. l2003j) 

2 Source terms 

In order to incorporate gravity we modified the original Relaxing TVD scheme 
through the addition of source terms within the Runge-Kutta algorithm. To 
achieve a good accuracy in simulations of near-hydrostatic equilibrium states, the 
gravity source terms are computed separately for the left- and right-going waves, 
by replacing zeros on the rhs. of eqns. (|1.5|) by the gravity source terms S(u L ) = 
{0,9xP L ,9yP L ,9ip L ,g L -™ L ) and S{u R ) = (0, g R p R , g R p R , g R p R , g R ■ m R ). The 
superscripts 'L' and 'R' in the gravitational acceleration reflects the fact that g 
is interpolated in a slightly different manner for the left-going and right -going 
waves. 



3 Constrained Transport 



The original RTVD MHD scheme by Pen et al. (120031) incorporates magnetic field 
evolution through the "constrained transport" (CT) algorithm (Evans et al. [1988]). 
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Magnetic field B is updated during the advection-constraint steps. The electro- 
motive force is computed in the advection step by the RTVD scheme and then 
used in the constraint step to preserve V • B = (Pen et al. 2003). 

Our extension of the RTVD algorithm includes an unsplit evolution of magneto- 
fluid, i.e. the simultaneous update of fluid variables and magnetic field in each 
Runge-Kutta step. 

4 Parallelization 

Piernik-MHD is fully parallelized with the aid of MPI library, by means of the 
block decomposition. Computational domain can be divided into any number of 
equal size blocks in any direction (see fig. [TJ. A test of code scalability has been 




Fig. 1. Total magnetic field in a global simulation of CR-driven galactic dynamo. The 
computational domain is dived into 1600 (20 x 20 x 4) equal blocks. 

done for HD Sedov explosion, in a series of experiments for different numbers of 
fixed size (n x — n y = n z = 64) MPI blocks, distributed over different CPU cores. 
The results displayed in Fig. [2] demonstrate that the growth of the core number 
from 8 up to 1024 results in the growth of the wall time only by a few percent. As 
it is apparent, the simplicity and homogeneity of the grid decomposition results 
in an excellent scaling of the code. 

5 Shearing box 

In addition to the the standard shearing model (Hawley et al. I1995|) . we imple- 
mented a modification of the shearing box method that allows to get rid of the 
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Piernik scaling (constant work per CPU) 
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Fig. 2. Results of the weak scaling test. In each test one CPU core processed one MPI 
block of fixed size (n x — n y = n z — 64). 



average shear velocity (v = —qfloz) when applying the CFL condition. Firstly, we 
make a simple transformation of initial conservative variables 

u = (p,pv x ,pv y ,pv z ,e) =$> u = (p,pv x ,pv y ,pv z ,e), (5.1) 

where v y = v y — v, e = e — + 2v y v) are quantities deprived of the terms 

containing mean shear flow. It can be shown that the transformation (|5.1|) does 
not break the conservative form of basic HD equations 

d t u + 8 x F{u, v x ) + d y G(u, v y ) + d z U(u, v z ) = S. (5.2) 

Following the fast Eulerian transport algorithm for differentially rotating disks 
introduced by Masset (MOO.) we can rewrite (|5.2[) as 

d t u + d x F(u, v x ) + d y G{u, v y ) + d z H(u, v z ) + vd y u = S. (5.3) 

The algorithm solving the equation (|5.3|) is then split into three steps: 

1. computation of source terms S(u), 

2. transport of u in x, y and z directions with v x ,v y , v z accordingly 

d t u + d x F(u, v x ) + d y G{u, v y ) + a 2 H(u, v.) = S, (5.4) 

using RTVD scheme, 

3. transport of u in y with v described by the linear advection equation 

d t u + vd y u — 0, (5.5) 
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The linear advection equation is being solved by the means of Fast Fourier Trans- 
form (FFT). Fluid variables are transformed into the frequency domain along y 
direction and phase shifted of <j> = vAt. The algorithm is now being extended to 
MHD equations in a manner preserving V • B = (Johnson et al. [2008) . 

6 Selfgravity 

We have implemented a Poisson solver in order to incorporate selfgravity of the 
fluids. Currently, our code supports selfgravity under condition of periodicity (or 
quasi-periodicity) of the domain in two or three directions, as our solver is based 
on FFT techniques. Utilisation of FFT is fully consistent with the previously 
described shearing box model. Poisson equation is solved in the shearing box by 

(1) the phase shift of the domain into the nearest periodic point (ID FFT in y), 

(2) application of suitable algorithm for periodic boundary conditions (two ID 
FFT in x), (3) shift back of the calculated gravitational potential (one ID FFT). 




Fig. 3. Snapshots of logarithm of surface density for times t = 3,5,9, 30fi _1 . Initially 
marginally stable state is slightly disturbed and undergoes fragmentation. Due to high 
cooling rate gas collapses to one gravitationally bounded object. 

The new, efficient shearing box algorithm has been tested in simulations of 



6 



Title : will be set by the publisher 



gravitational instability in protoplanetary disks fKowalik l2008j) . Following Gam- 
mie (|200ip we set initially uniform density distribution of ideal gas (7 = 2) with a 
subsonic (0.01c s ) velocity perturbation. Depending on how efficient the cooling of 
the gas is, fragmentation to gravitationally bounded object or gravitoturbulence 
may occur. Our results correspond closely to the results of a similar approach that 
uses a non-conservative scheme (Gammic 2001; Brandenburg & Dobler [2002p . 
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